ESTADISTICA PARA EL ANALISIS POLITICO II

Análisis de Conglomerados: Estrategia Basada en Densidad

Hasta ahora hemos encontrado clusters indicando cuantos se necesitaban, e indirectamente hemos forzado a que cada caso sea parte de uno de esos clusters. Veamos la data nuevamente:

# bibliotecas:
library(stringr)
library(magrittr)
library(htmltab)
library(factoextra)
library(cluster)

# coleccion
links=list(web="https://en.wikipedia.org/wiki/Democracy_Index",
           xpath ='//*[@id="mw-content-text"]/div/table[2]/tbody')
demo<- htmltab(doc = links$web, which =links$xpath)

# limpieza
names(demo)=str_split(names(demo),">>",simplify = T)[,1]%>%gsub('\\s','',.)
demo[,-c(1,8,9)]=lapply(demo[,-c(1,8,9)], trimws,whitespace = "[\\h\\v]")

# preparación
demo=demo[,-c(1)] #sin Rank
demo[,-c(1,8,9)]=lapply(demo[,-c(1,8,9)], as.numeric) # a numerico
row.names(demo)=demo$Country # cambiando row.names
demo=na.omit(demo)

# veamos que tenemos:
str(demo)
## 'data.frame':    167 obs. of  9 variables:
##  $ Country                        : chr  "Norway" "Iceland" "Sweden" "New Zealand" ...
##  $ Score                          : num  9.87 9.58 9.39 9.26 9.22 9.15 9.15 9.14 9.09 9.03 ...
##  $ Elec­toralpro­cessandplura­lism: num  10 10 9.58 10 10 9.58 9.58 10 10 9.58 ...
##  $ Functio­ningofgovern­ment      : num  9.64 9.29 9.64 9.29 9.29 7.86 9.64 8.93 8.93 9.29 ...
##  $ Poli­ticalpartici­pation       : num  10 8.89 8.33 8.89 8.33 8.33 7.78 8.33 7.78 7.78 ...
##  $ Poli­ticalculture              : num  10 10 10 8.13 9.38 10 8.75 8.75 8.75 9.38 ...
##  $ Civilliber­ties                : num  9.71 9.71 9.41 10 9.12 10 10 9.71 10 9.12 ...
##  $ Regimetype                     : chr  "Full democracy" "Full democracy" "Full democracy" "Full democracy" ...
##  $ Conti­nent                     : chr  "Europe" "Europe" "Europe" "Oceania" ...

Nuestro punto de partida clave siempre ha sido el cálculo de la matriz de distancias, añadamos la semilla aleatoria:

set.seed(2019)
inputData=demo[,c(3:7)]
g.dist = daisy(inputData, metric="gower")

Y con esa matriz calculamos cuatro clusters en nuestras clases previas, pero tal cantidad de clusters salió de una decisión algo arbitraria. Una pregunta exploratoria clave es cuántos clusters deberíamos calcular, y según ellos saber que hay una cantidad diferenciada de perfiles.

Determinando cantidad de clusters

Existe la medida gap, que sirve para determinar el mejor numero de clusters a pedir. Veamos:

fviz_nbclust(inputData, pam,diss=g.dist,method = "gap_stat",k.max = 10,verbose = F)

fviz_nbclust(inputData, hcut,diss=g.dist,method = "gap_stat",k.max = 10,verbose = F)

El numero de clusters varía, pero quedemonos con seis en ambos enfoques.

Evaluando los clusters obtenidos

Clusterizemos:

res.pam = pam(g.dist,6,cluster.only = F)
res.agnes = hcut(g.dist, k = 6,hc_func='agnes',hc_method = "ward.D")
res.diana = hcut(g.dist, k = 6,hc_func='diana')

Para evaluar, podemos analizar las siluetas (silhouettes), una medida que indica la calidad de asignación de un caso particular.

Evaluación gráfica

  • Evaluación gráfica para pam:
fviz_silhouette(res.pam)
##   cluster size ave.sil.width
## 1       1   14          0.56
## 2       2   29          0.25
## 3       3   31          0.34
## 4       4   35          0.09
## 5       5   32          0.26
## 6       6   26          0.28

  • Evaluación gráfica para agnes:
fviz_silhouette(res.agnes)
##   cluster size ave.sil.width
## 1       1   14          0.57
## 2       2   33          0.26
## 3       3   34          0.25
## 4       4   34          0.08
## 5       5   19          0.34
## 6       6   33          0.25

  • Evaluación gráfica para diana:
fviz_silhouette(res.diana)
##   cluster size ave.sil.width
## 1       1   39          0.47
## 2       2   52          0.27
## 3       3    3          0.25
## 4       4   38          0.10
## 5       5   18          0.34
## 6       6   17          0.25

Se asume que el gráfico que tiene menos siluetas negativas es el preferible a los demás.

Evaluación numérica

Esta etapa es para identificar a los casos mal asignados: los que tienen silueta negativa. Para ello es bueno saber lo que cada resultado nos trajo:

# por ejemplo tiene:
str(res.pam)
## List of 9
##  $ medoids   : chr [1:6] "Finland" "Cape Verde" "Guyana" "Benin" ...
##  $ id.med    : int [1:6] 8 26 55 82 126 155
##  $ clustering: Named int [1:167] 1 1 1 1 1 1 1 1 1 1 ...
##   ..- attr(*, "names")= chr [1:167] "Norway" "Iceland" "Sweden" "New Zealand" ...
##  $ objective : Named num [1:2] 0.0831 0.0794
##   ..- attr(*, "names")= chr [1:2] "build" "swap"
##  $ isolation : Factor w/ 3 levels "no","L","L*": 1 1 1 1 1 1
##   ..- attr(*, "names")= chr [1:6] "1" "2" "3" "4" ...
##  $ clusinfo  : num [1:6, 1:5] 14 29 31 35 32 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr [1:5] "size" "max_diss" "av_diss" "diameter" ...
##  $ silinfo   :List of 3
##   ..$ widths         : num [1:167, 1:3] 1 1 1 1 1 1 1 1 1 1 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:167] "Finland" "Iceland" "Denmark" "Sweden" ...
##   .. .. ..$ : chr [1:3] "cluster" "neighbor" "sil_width"
##   ..$ clus.avg.widths: num [1:6] 0.5571 0.2517 0.3405 0.0867 0.2638 ...
##   ..$ avg.width      : num 0.266
##  $ diss      : NULL
##  $ call      : language pam(x = g.dist, k = 6, cluster.only = F)
##  - attr(*, "class")= chr [1:2] "pam" "partition"

Aquí sabemos que res.pam es una lista con varios elementos. Uno de ellos es la información de siluetas, el cual tiene otros componentes:

str(res.pam$silinfo)
## List of 3
##  $ widths         : num [1:167, 1:3] 1 1 1 1 1 1 1 1 1 1 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:167] "Finland" "Iceland" "Denmark" "Sweden" ...
##   .. ..$ : chr [1:3] "cluster" "neighbor" "sil_width"
##  $ clus.avg.widths: num [1:6] 0.5571 0.2517 0.3405 0.0867 0.2638 ...
##  $ avg.width      : num 0.266

El primero, los widths, es donde esta la información de cada uno de los casos:

# veamos solo algunos.
head(res.pam$silinfo$widths)
##             cluster neighbor sil_width
## Finland           1        2 0.6829425
## Iceland           1        2 0.6818799
## Denmark           1        2 0.6763410
## Sweden            1        2 0.6637676
## Australia         1        2 0.6162141
## New Zealand       1        2 0.6122470

Creemos un data frame:

poorPAM=data.frame(res.pam$silinfo$widths)
poorPAM$country=row.names(poorPAM)

Nos interesa sólo sil_width negativos:

poorPAMcases=poorPAM[poorPAM$sil_width<0,'country']
# osea:
poorPAMcases
##  [1] "Senegal"          "Bhutan"           "Sri Lanka"       
##  [4] "Indonesia"        "Bangladesh"       "Papua New Guinea"
##  [7] "Turkey"           "Hong Kong"        "Guatemala"       
## [10] "Montenegro"       "Honduras"         "Namibia"         
## [13] "Belarus"          "Eswatini"         "Ethiopia"

La cantidad de paises es:

length(poorPAMcases)
## [1] 15

Podemos hacer lo mismo para las demás estrategias:

# agnes
poorAGNES=data.frame(res.agnes$silinfo$widths)
poorAGNES$country=row.names(poorAGNES)
poorAGNEScases=poorAGNES[poorAGNES$sil_width<0,'country']

#diana:
poorDIANA=data.frame(res.diana$silinfo$widths)
poorDIANA$country=row.names(poorDIANA)
poorDIANAcases=poorDIANA[poorDIANA$sil_width<0,'country']

Ahora que tenemos todos los paises mal asignados, podemos interrogar a los resultados usando teoría de conjuntos, por ejemplo:

  • Los paises mal asignados en agnes y en pam:
intersect(poorAGNEScases,poorPAMcases)
## [1] "Hong Kong"  "Bhutan"     "Bangladesh" "Honduras"   "Guatemala"
  • Los paises mal asignados por agnes pero no por pam:
setdiff(poorAGNEScases,poorPAMcases)
##  [1] "Uruguay"      "Mauritius"    "South Africa" "Malaysia"    
##  [5] "Benin"        "Georgia"      "Ivory Coast"  "Fiji"        
##  [9] "Myanmar"      "Niger"        "Jordan"       "Kuwait"      
## [13] "Zimbabwe"     "Angola"
  • Los paises mal asignados por pam pero no por agnes:
setdiff(poorPAMcases,poorAGNEScases)
##  [1] "Senegal"          "Sri Lanka"        "Indonesia"       
##  [4] "Papua New Guinea" "Turkey"           "Montenegro"      
##  [7] "Namibia"          "Belarus"          "Eswatini"        
## [10] "Ethiopia"
  • Los paises mal asignados por pam o por agnes:
union(poorPAMcases,poorAGNEScases)
##  [1] "Senegal"          "Bhutan"           "Sri Lanka"       
##  [4] "Indonesia"        "Bangladesh"       "Papua New Guinea"
##  [7] "Turkey"           "Hong Kong"        "Guatemala"       
## [10] "Montenegro"       "Honduras"         "Namibia"         
## [13] "Belarus"          "Eswatini"         "Ethiopia"        
## [16] "Uruguay"          "Mauritius"        "South Africa"    
## [19] "Malaysia"         "Benin"            "Georgia"         
## [22] "Ivory Coast"      "Fiji"             "Myanmar"         
## [25] "Niger"            "Jordan"           "Kuwait"          
## [28] "Zimbabwe"         "Angola"

Density Based clustering

La estrategia basada en densidad sigue una estrategia muy sencilla: juntar a los casos cuya cercanía entre sí los diferencia de otros. Aquí hay un resumen breve del tema:

El algoritmo dbscan requiere dos parametros:

  1. La distancia epsilon a usar para clusterizar los casos.
  2. La cantidad k minima de puntos para formar un cluster.

El valor k que se usará es al menos la cantidad de dimensiones. Como se han usado cinco variables, usaremos k=5. Calculemos el epsilon.

Sin embargo, el principal problema es que necesitamos un mapa de posiciones para todos los paises. Eso requiere una técnica que proyecte las dimensiones originales en un plano bidimensional. Para ello usaremos la técnica llamada escalamiento multidimensional:

proyeccion = cmdscale(g.dist, k=2,add = T) # k is the number of dim

Habiendo calculado la proyeccción, recuperemos las coordenadas del mapa del mundo basado en nuestras dimensiones nuevas:

# data frame prep:
inputData$dim1 <- proyeccion$points[,1]
inputData$dim2 <- proyeccion$points[,2]

Aquí puedes ver el mapa:

base= ggplot(inputData,aes(x=dim1, y=dim2,label=row.names(inputData))) 
base + geom_text(size=2)

Obtengamos los clusters anteriores:

inputData$pam=as.factor(res.pam$clustering)
inputData$agnes=as.factor(res.agnes$cluster)
inputData$diana=as.factor(res.diana$cluster)

Grafiquemos:

# Estimado limites:
min(inputData[,c('dim1','dim2')]); max(inputData[,c('dim1','dim2')])
## [1] -0.6287386
## [1] 0.6412899
limites=c(-0.7,0.7)
base= ggplot(inputData,aes(x=dim1, y=dim2)) + ylim(limites) + xlim(limites) + coord_fixed()
base + geom_point(size=2, aes(color=pam))  + labs(title = "PAM") 

base + geom_point(size=2, aes(color=agnes)) + labs(title = "AGNES")

base + geom_point(size=2, aes(color=diana)) + labs(title = "DIANA")

Ahora calculemos usando dbscan:

  1. Nuevas distancias
# euclidea!!
g.dist.cmd = daisy(inputData[,c('dim1','dim2')], metric = 'euclidean')
  1. Calculo de epsilon
library(dbscan)
kNNdistplot(g.dist.cmd, k=5)

  1. Obteniendo clusters
library(fpc)
db.cmd = dbscan(g.dist.cmd, eps=0.06, MinPts=5,method = 'dist');  db.cmd
## dbscan Pts=167 MinPts=5 eps=0.06
##         0  1  2  3
## border 11  2  0 12
## seed    0 49 33 60
## total  11 51 33 72

Aqui asignamos los clusters a otra columna:

inputData$dbCMD=as.factor(db.cmd$cluster)
  1. Graficando

Aquí sin texto:

library(ggrepel)
base= ggplot(inputData,aes(x=dim1, y=dim2)) + ylim(limites) + xlim(limites) + coord_fixed()
dbplot= base + geom_point(aes(color=dbCMD)) 
dbplot

Aquí con mucho texto:

dbplot + geom_text_repel(size=5,aes(label=row.names(inputData)))

Aquí sólo los atípicos:

LABEL=ifelse(inputData$dbCMD==0,row.names(inputData),"")
dbplot + geom_text_repel(aes(label=LABEL),
                         size=5, 
                         direction = "y", ylim = 0.45,
                         angle=45,
                         segment.colour = "grey")

Nota que en esta técnica hay casos que no serán clusterizados.


al INICIO

VOLVER A CONTENIDOS